Response of discrete nonlinear systems with many degrees of freedom 
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We study the response of a large array of coupled nonlinear oscillators to parametric excitation, 
motivated by the growing interest in the nonlinear dynamics of microelectromechanical and nano- 
electromechanical systems (MEMS and NEMS). Using a multiscale analysis, we derive an amplitude 
equation that captures the slow dynamics of the coupled oscillators just above the onset of para- 
metric oscillations. The amplitude equation that we derive here from first principles exhibits a 
wavenumber dependent bifurcation similar in character to the behavior known to exist in fluids un- 
dergoing the Faraday wave instability. We confirm this behavior numerically and make suggestions 
for testing it experimentally with MEMS and NEMS resonators. 
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I. MOTIVATION 

In the last decade we have witnessed exciting tech- 
nological advances in the fabrication and control of mi- 
croelectromechanical and nanoelectromechanical systems 
(MEMS and NEMS). Such systems are being developed 
for a host of nanotechnological applications, as well as 
for basic research in the mesoscopic physics of phonons, 
and the general study of the behavior of mechanical de- 
grees of freedom at the interface between the quantum 
and the classical worlds 

[3, HQ. Surprisingly, NEMS 
have also opened up a new experimental window into 
the study of the nonlinear dynamics of discrete systems 
with many degrees of freedom. A combination of three 
properties of NEMS resonators has led to this unique ex- 
perimental opportunity. First and most important is the 
experimental observation that micro- and nanomechani- 
cal resonators tend to behave nonlinearly at very modest 
amplitudes. This nonlinear behavior has not only been 
observed experimentally 0, 0, IE S IE S E3, IlJ. but 
has already been exploited to achieve mechanical signal 
amplification and mechanical noise squeezing [f3 . Il3| in 
single resonators. Second is the fact that at their di- 
mensions, the normal frequencies of nanomechanical res- 
onators are extremely high — recently exceeding the 1GHz 
mark 0, — facilitating the design of ultra-fast me- 
chanical devices, and making the waiting times for un- 
wanted transients bearable on experimental time scales. 
Third is the technological ability to fabricate large ar- 
rays of MEMS and NEMS resonators whose collective 
response might be useful for signal enhancement and 
noise reduction ^(|, as well as for sophisticated me- 
chanical signal processing applications. Such arrays have 
already exhibited interesting nonlinear dynamics rang- 
ing from the formation of extended patterns — as 
one commonly observes in analogous continuous systems 



such as Faraday waves — to that of intrinsically localized 
modes 0,E1|2(J. Thus, nanomechanical resonator ar- 
rays are perfect for testing dynamical theories of discrete 
nonlinear systems with many degrees of freedom. At the 
same time, the theoretical understanding of such systems 
may prove useful for future nanotechnological applica- 
tions. 
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Two of us (Lifshitz and Cross 0, henceforth LC]) have 
recently studied the response of coupled nonlinear oscil- 
lators to parametric excitation. We used secular per- 
turbation theory to convert the equations of motion for 
the oscillators into a set of coupled nonlinear algebraic 
equations for the normal mode amplitudes of the sys- 
tem, enabling us to obtain exact results for small arrays 
but only a qualitative understanding of the dynamics of 
large arrays. In order to obtain analytical results for large 
arrays we study here the same system of equations, ap- 
proaching it from the continuous limit of infinitely-many 
degrees of freedom. A novel feature of the parameri- 
cally driven instability is that the bifurcation to stand- 
ing waves switches from supercritical (second order) to 
subcritical (first order) at a wave number at or close to 
the critical one for which the required driving force is 
minimum. This changes the form of the amplitude equa- 
tion that describes the onset of the paramerically driven 
waves so that it no longer has the standard "Ginzburg- 
Landau" form. Our central result is this new scaled am- 
plitude equation (|12fl . governed by a single control pa- 
rameter, that captures the slow dynamics of the coupled 
oscillators just above the onset of parametric oscillations 
including this unusual bifurcation behavior. We confirm 
the behavior numerically and make suggestions for test- 
ing it experimentally. Although our focus is on para- 
metrically driven oscillators, the amplitude equation we 
derive should also apply to other parametrically driven 
wave systems with weak nonlinear damping. 
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II. EQUATIONS OF MOTION 

The equations of motion introduced in LC are 

3 1 , 

u„ + u„ + u n - -r(u n+ i - 2u n + u„_i) 

+ ^A 2 [l + H COs(2uj p t)] (Un+l - 2u n + ttn-l) 

- ^v[( u n+l - U n ) 2 (u n+1 ~ il n ) 

- {u n - u„_i) 2 (w„ - u„_i)] = 0, (1) 

with n = 1 . . . N, and fixed boundary conditions uo = 
Ujv+i = 0. The equations of motion are modeled after 
the experiment of Buks and Roukes [17| , who succeeded 
in fabricating, exciting, and measuring the response to 
parametric excitation of an array of 67 micromechani- 
cal resonating gold beams. Detailed arguments for the 
choice of terms introduced into the equations of motion 
are given by LC. The guiding principle is to introduce 
only those terms that are essential for capturing the phys- 
ical behavior observed in the experiment. These include a 
cubic nonlinear elastic restoring force (whose coefficient 
is scaled to 1), a dc electrostatic nearest-neighbor cou- 
pling term with a small ac component responsible for 
the parametric excitation (with coefficients A 2 and A 2 H 
respectively), and linear as well as cubic nonlinear dis- 
sipation terms (with coefficients T and 77 respectively). 
The nonlinear contribution to the dissipation, though ex- 
pected to be small, is essential for the saturation of the 
amplitude of the parametrically driven motion |2lJ and is 
important in the behavior of the bifircation as a function 
of wave number. Both dissipation terms are taken to be 
of a nearest neighbor form, motivated by the experimen- 
tal indication that most of the dissipation comes from the 
electrostatic interaction between neighboring beams. It 
should be noted that the effect of introducing gradient- 
dependent dissipation terms instead of local dissipation 
terms is merely to renormalize the bare dissipation coef- 
ficients r and 77 to wavenumber-dependent coefficients of 
the form Tsui 2 (q/2) and 77 sin 2 (g/2) respectively, as will 
become apparent below. 

The dissipation of the system is assumed to be weak, 



which makes it possible to excite the beams with rel- 
atively small driving amplitudes. In such case the re- 
sponse of the beams is moderate, justifying the descrip- 
tion of the system with nonlinearities up to cubic terms 
only. The weak dissipation can be parameterized by in- 
troducing a small expansion parameter e <C 1, physically 
defined by the linear dissipation coefficient T = £7, with 
7 of order one. The driving amplitude is then expressed 
by A 2 H = eh, with h of order one. We assume the 
system is excited in its first instability tongue, i.e. we 
take Up itself to lie within the normal frequency band 
y/1 — 2 A 2 < ui p < 1. The weakly nonlinear regime is 
studied by expanding the displacements u n in powers of 
e. Taking the leading term to be of the order of e 1 / 2 
ensures that all the corrections, to a simple set of equa- 
tions describing N coupled harmonic oscillators, enter 
the equations at the same order of e 3 / 2 . 



III. AMPLITUDE EQUATIONS FOR COUNTER 
PROPAGATING WAVES 

We introduce a continuous displacement field u(x,t), 
keeping in mind that only for integral values x = n of 
the spatial coordinate does it actually correspond to the 
displacements u(n, t) = u n {t) of the discrete set of oscilla- 
tors in the array. We introduce slow spatial and temporal 
scales, X = ex and T = et, upon which the dynamics of 
the envelope function occurs, and expand the displace- 
ment field in terms of e, 

u(x,t) = e 1/2 [(A + {X,T)e~ iq ? x + A*_(X,T)e iq * x ) e 1 ""* 
+ c.c] + e 3 / 2 7i« (x, t, X, T) + . . . , (2) 

where the asterisk and c.c. stand for the complex conju- 
gate, and q p and u p are related through the dispersion re- 
lation uj 2 — 1— 2A 2 sin 2 (<7 p /2). The response to lowest or- 
der in e is expressed in terms of two counter-propagating 
waves with complex amplitudes A + and A_, a typical 
ansatz for parametrically excited continuous systems |22| . 
We substitute the ansatz J3J into the equations of mo- 
tion (JIJ term by term. Up to order e 3 / 2 , we have 



i 1/2 lu 2 v (A+ (X, T)e- iq " x + A*_ {X, T)e lq ? x ) e*"'* + e 3/2 ^^(x, t, X, T) 

3/2o- fdA+(X,T) _ ia x dA*_(X,T) . \ lu t 
+ e i/2 2tuj v [ — -e lQpX H ^ — -e lq " x ) e luJpt + c.c. ; 



dT dT 
u n±l = e 1 ' 2 {A + {X,T)e- iq ^ x±v > + A*_ (X,T)e iq ^ x±1 A e iuj ^ + e^ 2 u w (x±l,t,X,T) 

1 3/2 ( 9A+ -iq p x p =F»gy _1_ iq v x ±iq„\ iuj p t , . 

\dX + dX ) + ' 

^A 2 K +1 - 2u n + u n -x) - -e 1 / 2 2A 2 sin 2 (g p /2) {A + e~ iq " x + A*_e iq » x ) 



(3a) 



(3b) 
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+ e 3 ' 2 (u^ [x + l,t, X, T) — 2vfV {x, t, X, T) + u^(x -l,t,X, T)) 



e 3 / 2 A 2 isin(g p ) (^±e 



■iq p x 



dA*_ 
~dX 



c. ; 



lehcos{2cj p t)(u n+1 - 2u n + w„_ x ) = -e 3 / 2 /i.sin 2 (<7 p /2) (A-e^** + A* + e lqpX ) e* 1 ** + 0(e l3uJpt ) + 



c.c. ; 



;e 7 (it„ + i - 2u„ + u n _i) = -e 3 / 2 2i 7 w p sin 2 (g p /2) (A+e"^ + Ale* 9 "*) e*"»* + c.c. ; 

3 =e 3/2 3 [(|A+| 2 + 2\A^\ 2 )A +e - iq " x + {2\A+\ 2 + \A^\ 2 )A*_e iq r x ] e Wpt + O (e^'.e* 3 * 1 ) + c.c. : 



[ ( u n+ i - u n ) 2 (ii n+ i - u n ) - (u n - u n -i) 2 {u n - u n - 



i8r]ujp sin 4 (q p /2) 



x [(\A+\ 2 + 2\A^\ 2 )A+e- tq " x + (2\A + \ 2 + \A_\ 2 )A*_e iq » x ] e*"'* + 0{e l3u)p \ e l3qpX ) + c.c. 

I 



(3c) 
(3d) 

(3e) 
(3f) 

(3g) 



where 0{e l3u)pt , e l3q " x ) are terms proportional to e l3ulpt 
or e l3qpX which do not enter the dynamics at the low- 
est order of the e expansion. At the order of e 1 / 2 , the 
equations of motion Q are satisfied trivially, yielding 
the dispersion relation mentioned earlier. At the order 
of e 3 / 2 on the other hand we must apply a solvability 
condition |22J , which requires that all terms proportional 
to e^pti^pz) mus t vanish. We therefore obtain the two 
coupled amplitude equations, 



dA± dA± 

± v n 



dT 



dX 



-7 sin 



m — A± =F * sin — 

V 2 / 2u p V 2 / 



A 3 



4??Sin4 (l) Ti 4 



I Ah 



2\A T \ 2 ) A ± ,(4) 



where the upper signs (lower signs) give the equation 
for A+ (A_) obtained from the restriction on the terms 
proportional to e^ UJpt ~ qpX ' ) ( e ( lw p*+9p a: )) j an d 



dq 



A 2 sin(tfr) 
2uj v 



(5) 



is the group velocity. A detailed derivation of the am- 
plitude equations Q can be found in Ref. I23L Similar 
equations were previously derived for describing Faraday 
waves I2II2H. 



IV. 



REDUCTION TO A SINGLE AMPLITUDE 
EQUATION 



By linearizing the amplitude equations (J2J about the 
zero solution (A + = A_ = 0) we find that the linear 
combination of the two amplitudes that first becomes 
unstable at h — 2 r yuj p is B oc (A_|_ — iAJ) — representing 
the emergence of a standing wave with a temporal phase 
of 7r/4 relative to the drive — while the orthogonal linear 
combination of the amplitudes decays exponentially and 
does not participate in the dynamics at onset. Thus, 
just above threshold we can reduce the description of the 
dynamics to a single amplitude B, where at a finite am- 
plitude above threshold a band of unstable modes around 
q p can contribute to the spatial form of B. This is similar 



to the procedure introduced by Riecke [26| for describing 
the onset of Faraday waves. 

In the absence of nonlinear damping 7/ = the coeffi- 
cient of the nonlinear term in Eq. (0J is purely imaginary. 
This turns out to yield the result that there is no non- 
linear saturation term of the usual form \B\ 2 B for the 
standing waves at the resonant wave vector q p . An anal- 
ysis of the saturation of standing waves of nearby wave 
numbers q p + k for this case shows that the coefficient of 
\Bk\ 2 Bk is positive for k < 0, zero at k = 0, and negative 
for k > O.This means that for 77 = the nature of the bi- 
furcation to standing waves switches from supercritical to 
subcritical precisely at the critical wave number q = q p . 
Thus the amplitude equation describing the onset of the 
standing waves will not be of the usual Ginzburg-Landau 
form. It is the goal of the present section to derive the 
novel form of the amplitude equation. We will do this 
for the case of 77 small but non-zero, in which case the 
switch from subcritical to supercritical occurs at a wave 
number close to the resonant wave number q p . 

We define a reduced driving amplitude g with respect 
to the threshold 27^ by letting (h — 2 r yLu p )/2ju!p = gS, 
with 5 « 1. In order to obtain an equation, describ- 
ing the relevant slow dynamics of the new amplitude B, 
we need to select the proper scaling of the original ampli- 
tudes A±, as well as their spatial and temporal variables, 
with respect to the new small parameter 6. We assume 
that the coefficient of nonlinear dissipation 77 is small. As 
we have seen, for some wave numbers near q p the bifurca- 
tion is then subcritical, so that a quintic term must enter 
in order to saturate the growth of the amplitude. This 
is similar to the situation encountered by Deissler and 
Brand 27] who studied localized modes near a subcritical 
bifurcation to traveling waves. Here this can be achieved 
by defining the small parameter S with respect to the 
coefficient of nonlinear dissipation — letting 77 = S 1 / 2 ^, 
with r/Q of order one — and taking the amplitudes to be of 
order S 1 ^ 4 . Further noting that with a drive amplitude 
that scales as 5 the growth rate scales like 5 as well, and 
the bandwidth of unstable wave numbers scales as S 1 / 2 , 
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we finally make the ansatz that 



where O is a linear operator given by the matrix 



A, 
A. 



(6) 



5/4 / w W(X,T,£,t) 
~ " 1 vW(X,T,ir)J> 



where f = S 1 / 2 X and f = ST are the new spatial and 
temporal scales respectively. The amplitude equation for 
B(£,t) is derived by once again using the multiple scales 
method. We substitute the ansatz © into the coupled 
amplitude equations Q and collect terms of different or- 
ders of 8. 

Again, to the lowest order of expansion the equations 
are satisfied trivially. Collecting all terms of order J 3 / 4 
in Eqs. J3J yields 



dB 



+ t 



9 

2u> Tl 



\B\ 2 B 



1 

—4 



(7) 



dr + Vgd x + 7sin 2 (q p /2) i7sin 2 (<7 p /2) 

-i7sin 2 ((jp/2) <9 T - w 3 9 x + 7sin 2 (g p /2) 

There is no solvability condition to satisfy here because 
the right hand side of Eq. (JJJ is automatically orthogonal 

to the zero eigenvalue of O, ( j V The solution of Eq. Q 

is given by 



V W 

1 



dB 9 



2 7 sin 2 (gp/2) V <9£ 2w p 



(9) 



We plug Eq. © back into Eqs. (0J, collect all the terms 
of order 5 5 / 4 and obtain 







w (2) 
w (2) 



as 



a 2 s 



2 7 sin 2 (g p /2) d£ 2 



7sin 2 (g p /2) g S - 12^ sin 4 (g p /2)|B| 2 B 



+ i 



3(-« fl ) 



2w P 7sin (g p /2) 



4151 



.as 



81 



8w 2 7sin% p /2) 



|B| 4 S 



(10) 



The solvability condition of Eq. , stating that the zero 

eigenvalue of the operator D, ( * ^ , must be orthogonal 

to the right hand side of the equation, determines the 
amplitude equation for B |22| . Thus, the expression in 
square brackets in Eq. i|l(J|) must vanish. After applying 
one last set of rescaling transformations, 

f = j |_3 \v g \t 

32 t^r/^ s i n 10 (g p /2) ' 8 w p r?o7 sin 6 (g p /2) ' 

|S| 2 -^07Bin 6 (|)|B| 2 , 
^|^ 2 sin 8 (|) ff , (11) 

we end up with an amplitude equation governed by a 
single parameter, 

8B 8 2 B 2f A . nl2 dB n2 dB*\ 

-2\B\ 2 B-\B\ 4 B. (12) 

Note that Eq. (|12|l is a scaled equation with all coeffi- 
cients of order unity. In an equation for the unsealed 
mode amplitude B = 8 1 / i B (i.e. the combination of 
A + , A- that goes unstable) and using the length and 
time variables X, T rather than £, r, the coefficient of the 



cubic saturation term \B\ 2 B would be small, scaling as 
(5 1 / 2 , reflecting its origin in the nonlinear damping term, 
the linear growth term proportional to B would have a 
coefficient scaling as (h—h c )/h c , and all other coefficients 
would be of order unity. All terms in the equation are 
then of comparable size for (h—h c )/h c ~ S, for oscillation 
amplitudes B ~ 8 1 / 4 , and for length scales X <~ <5~ 4 / 2 
and time scales T ~ <5~ 4 . The form of Eq. 1121) is also 
applicable to the onset of parametrically driven standing 
waves in continuum systems with weak nonlinear damp- 
ing, and combines in a single e quation a number of effects 
studied previously [24L l25l 12a . I27L I28L l29j . The numerical 
coefficients in the variable scalings and in the equation 
will, of course, be different for different systems. 

Once we have obtained this amplitude equation it can 
be used to study a variety of dynamical solutions, ranging 
from simple single-mode to more complicated nonlinear 
extended solutions, or possibly even localized solutions. 



V. SINGLE MODE OSCILLATIONS 

In this section we begin the investigation of the so- 
lutions of Eq. (|12|l by focusing on the regime of small 
reduced drive amplitude g and look upon the saturation 
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FIG. 1: (Color online) Response of the oscillator array plot- 
ted as a function of reduced amplitude g for three different 
scaled wave number shifts determined by fixing the number 
of oscillators N: k = (for N = 100) and k = -0.81 (for 
N — 92), which bifurcate supercritically, and k — 1.55 (for 
N = 98) which bifurcates subcritically showing clear hys- 
teresis. Solid and dashed lines are the positive and negative 
square root branches of the calculated response in (1141 . the 
latter clearly unstable. Open circles are numerical values ob- 
tained by integration of the equations of motion Q for each 
N, with A = 0.5, Lu p = 0.767445, e-y = 0.01, and t) = 0.1. 
For these parameters a unity of the scaled squared amplitude 
|fefe| 2 corresponds to 2 10 -4 for the unsealed squared ampli- 
tude, and a unity of the wave number shift k corresponds to 
a shift of 9 10 -3 of the unsealed wavenumber. Dots show the 
subcritically-bifurcating single- mode solution of LC [2l| . 

of single-mode solutions of the form 

B = b k e-^, (13) 

corresponding — via the scaling £ = S%£, where the scale 
factor 5j is defined in l|ll|) — to a standing wave with a 
shifted wave number q = q p + ke^/S/S^. From the lin- 
ear terms in the amplitude equation (|12f> we find, as ex- 
pected, that for g > k 2 the zero-displacement solution 
is unstable to small perturbations of the form of ljl3|l . 
defining the parabolic neutral stability curve, shown as 
a dashed line in Fig. [5] The nonlinear gradients and 
the cubic term take the simple form 2(k — l)\bk\ 2 bk- For 
k < 1 these terms immediately act to saturate the growth 
of the amplitude assisted by the quintic term. Standing 
waves therefore bifurcate supercritically from the zero- 
displacement state. For k > 1 the cubic terms act to 
increase the growth of the amplitude, and saturation is 
achieved only by the quintic term. Standing waves there- 
fore bifurcate subcritically from the zero-displacement 
state. Note that wave-number dependent bifurcations 
similar in character were also predicted and observed nu- 
merically in Faraday waves [25L l28l |29| . The saturated 
amplitude \bk\, obtained by setting Eq. ljl"2|) to zero, is 
given by 

\b k \ 2 = (k - 1) ± V(* - I) 2 + (9 - k 2 ) > 0. (14) 



The original boundary conditions u(0, t) = u(N + 1, t) = 
impose a phase of 7r/4 on bk and require that the wave 
numbers be quantized q m — rmr/(N + 1), m = 1 . . . N. 

In Fig. (JIJ we plot \bk\ 2 as a function of the reduced 
driving amplitude g for three different wave number shifts 
k. The solid (dashed) lines are the stable (unstable) so- 
lutions of Eq. I|14|) . The circles were obtained by numer- 
ical integration of the equations of motion QJ. For each 
driving amplitude, the Fourier components of the steady 
state solution were computed to verify that only single 
modes are found, suggesting that in this regime of pa- 
rameters only these states are stable. By changing the 
number N of oscillators we could control the wave num- 
ber shift k for a fixed value of u> p . In experiment it might 
be easier to control k, for a fixed value of u p , by changing 
the dc component of the potential difference between the 
beams, thus changing the dispersion relation and with it 
the value of q p . 

Lifshitz and Cross 0, Eq. (33)] obtained the ex- 
act form of single mode solutions by substituting u n — 
A m sm{q m n)e lu ' pt + c.c directly into the equations of 
motion Q). In the limit of driving amplitudes just 
above threshold and 77 <C 1 their solution corresponds 
to Eq. I|14|) . as shown by the dots in Fig. QJ. In or- 
der to compare the two solutions one should note that 
in both cases the system oscillates in one of its normal 
modes q m with the driving frequency cu p . Here we use k 
to denote the difference between q p and the wavenumber 
of the oscillating pattern, whereas Lifshitz and Cross use 
a frequency detuning oj p — ui m to denote the difference 
between the normal frequency and the actual frequency 
of the oscillations. The frequency detuning lu p — uj m is 
proportional to v g ke^/S, implying that for an infinitly 
extended system the standing waves will always bifur- 
cate supercritically with a wave number q p if the driving 
amplitude is increased quasistatically. It is the discrete- 
ness of the normal modes which provides the detuning 
essential for a subcritical bifurcation if only quasistatic 
changes are performed. 



VI. SECONDARY INSTABILITIES 

We study secondary instabilities of the single mode 
solutions by performing linear stability analysis of 11311 . 
We substitute 

B(£, T ) = 6 fe e- 4fc « + (/3 + (r)e- l(fc+Q) « +^I(r)e- j(fc - Q) «) , 

(15) 

with \(3±\ <C 1, into the amplitude equation 1121) and 
linearize in j3±. Since the amplitude equation (|12H is 
invariant under phase transformations B —> Be~ ltp , the 
stability of the single mode solution cannot depend on 
the phase of bk- We therefore assume that bk is real, and 
linearize the following terms of Eq. (|12fl 
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\B\*B - e" lfc « (&3 +b 2 k ((2P + +p_)e- i( X + (/?; + 2/3* )e^)) ; 
|S| 4 B -► e~ lfc « (bf + b 4 k ((3/3+ + 2/3_)e^ Q « + (2/?; + 3/3* )e lQ «)) ; 

S 2 — -> e- lfc « (^ + i62((2A:/3 + + (fc-Q)/3_)e^ + ((fc + Q)/3; + 2fc/3*)e l< 3«)). 



(16) 



The terms of order one of the equation obtained from the linearization of Eq. I|12fl recover the same Eq. i|14fl for the 
single mode amplitude b k . The terms with spatial dependence of e~ 1 ^ and e 1 ^ determine the temporal development 
of the perturbations /3+ and /3_ respectively, and are given by 

3 ^-m(M, ,17) 



M 



2b 2 (k - 1 - b 2 k ) -Q 2 - 2Q{k - 46 2 /3) 2b\{k -l-b 2 + Q/3) 

2b 2 (k -l-b 2 - Q/3) 2b 2 (k -l-b 2 )-Q 2 + 2Q(k - 46 2 /3) 



where we have used Eq. (|14|l 

g _ fc 2 + 4(fc - 1)6 2 - 3b 4 k = 2b\(k - 1 - b\), (18) 

The single mode solution (|13|) is a stable solution only 
if all wave numbers Q yield negative eigenvalues for M, 
obtained for 

trM = 46|(fc — 1 — b 2 k ) - 2Q 2 < and (19) 
\M\ = 8 -Q 2 Qq 2 bi + 5 -^b 2 h 2 } > 

for all wave numbers Q. The negative square root branch 
in i|14|) (obtained for b 2 < k — 1) is confirmed to be 
always unstable. The stability of the positive square 
root branch, determined by setting the inequalities in 
Eqs. (|19fl to equalities, is bounded by the dotted curve 
in Fig. |21 describing the stability balloon of the single 
mode state. Outside the stability balloon the standing 
wave undergoes an Eckhaus instability with wave num- 
bers k ± Q, which occurs first at Q — > 0. When taking 
into consideration the discreteness of the system, where 
only wave numbers with Q > AQjy — S^n / e^/~5{N + 1) 
can be taken into account, the stability balloon is ex- 
tended to the solid curve in the Figure, plotted for the 
case of N = 92. 

A transition from one single-mode state to a new 
single- mode state with a wave number shift of nAQjv, 
for some integer n, occurs once the driving amplitude is 
increased and has crossed the upper bound of the sta- 
bility balloon. Since the upper bound monotonically 
increases with k, the new wave number will always be 
larger. A sequence of three transitions, obtained numer- 
ically, is shown in Fig. [21 superimposed with our theo- 
retical predictions. The sequence of transitions is also 
sketched for comparison within the stability balloon in 
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FIG. 2: (Color online) Stability boundaries of the single- 
mode solution of Eq. 112H in the g vs. k plane. Dashed line: 
neutral stability boundary. Dotted line: stability boundary 
of the single-mode solution 1131 for a continuous spectrum. 
Solid line: stability boundary of the single-mode solution for 
N — 92 and the parameters of Fig. Q For k > 1 the bi- 
furcation from zero displacement becomes subcritical and the 
lower stability boundary is the locus of saddle-node bifurca- 
tions (green line). Vertical and horizontal arrows mark the 
secondary instability transitions shown in Fig. |3] 



Fig. This type of analysis yields predictions for hys- 
teretic transitions on slow sweeps, with results for when 
the transition occurs and for the new state that devel- 
ops, which for larger g must be selected out of a band 
of possible stable states. It is clear that secondary in- 
stabilities allow an additional scenario for observing sub- 
critical transitions, even in very large systems, where the 
wavenumber separation is small. Once a secondary tran- 
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FIG. 3: (Color online) A sequence of secondary instabilities 
following the initial onset of single-mode oscillations in an 
array of 92 beams, with the parameters of Fig.0 plotted as a 
function of the reduced driving amplitude g. Solid (dashed) 
lines are stable (unstable) solutions defined by 1141 . for k — 
—0.88, k' = 2.81, and k" — 6.51, corresponding to the first 
wave number k to emerge and two shifts to k + AQ m and k + 
2AQm respectively. Numerical integration of the equations 
of motion Q for an upward sweep of g (blue x's), followed 
by a downward sweep (red o's) exhibits clear hysteresis and 
confirms the theoretical predictions for the stability of single- 
mode oscillations as illustrated in Fig. 01 

sition has occurred, the system will return to oscillate in 
the original wavenumber only when reducing the driving 
amplitude below the saddle node point of the subcritical 
bifurcation. 



VII. CONCLUSIONS 

We derived amplitude equations (@J describing the re- 
sponse of large arrays of nonlinear coupled oscillators 
to parametric excitation, directly from the equations of 
motion yielding exact expressions for all the coefficients. 
The dynamics at the onset of oscillations was studied by 
reducing these two coupled equations into a single scaled 
equation governed by a single control parameter (|12|l . 
Single mode standing waves were found to be the initial 
states that develop just above threshold, typical of para- 



metric excitation. The single mode oscillations bifurcate 
from the zero-displacement state either supercritically or 
subcritically, depending on the wave number of the os- 
cillations. The wave number dependence originates from 
the nonlinear gradient terms of the amplitude equation, 
which were usually disregarded in the past by others in 
the analysis of parametric oscillations above threshold. 
We also examined the stability of single mode oscilla- 
tions, predicting a transition of the initial standing wave 
state to a new standing wave with a larger wave number 
as the driving amplitude is increased. 

In this work we showed that interesting response of 
coupled nonlinear oscillators excited parametrically can 
also be obtained for quasistatic driving amplitude sweeps, 
rather than the frequency sweeps that are usually pre- 
formed in experiments. We proposed, and numerically 
demonstrated, two experimental schemes for observing 
our predictions, hoping to draw more experimental at- 
tention to the dynamics produced by such sweeps of the 
driving amplitude. 

The results obtained by the numerical integration of 
the equations of motion agree with our analysis, sup- 
porting the validity of the amplitude equation (|12f) . We 
therefore believe that the amplitude equations we de- 
rived can serve as a good starting point for studying 
other possible states of the system. One particular in- 
teresting dynamical behavior that can be considered is 
that of localized modes, often observed in arrays of cou- 
pled nonlinear oscillators and in other nonlinear systems 
as well. The conditions for obtaining such modes and 
their dynamical properties could be studied by looking 
for localized states of the amplitude equations. Another 
interesting aspect that can be addressed using the ampli- 
tude equations we derived is the response of the array to 
fast (rather than quasistatic) sweeps of the driving am- 
plitude, which should lead to more complicated nonlinear 
response as observed by LC in their work. 
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